class IWLS_G
{
	public:

	IWLS_G(void)
	{	
		boost::math::normal_distribution<> *d=new boost::math::normal_distribution<>(0,1); 
		_d=d;
	}
	~IWLS_G()
	{
		delete _d;
	}
	void operator()(mat X,double *Y,mat m, double si,int nu=1)
	{
		double e=1;
		double thres=0.0001;
		int n=X.n_rows;
		int p=X.n_cols;
		mat s(p,p);
	        s.fill(si);
		s=diagmat(s);	
		mat Yt=vect2mat<double>(Y,n,1);
		mat I;
		I.eye(p,p);
		mat Xt=rbind<mat>(X,I);
		mat b=inv(X.t()*X)*X.t()*Yt;
		mat V(p,p);
		mat z(n,1);
		mat w(n,1);
		int i=0;
		while(e>thres)
		{
			mat pn=pnorm(X*b);
			mat dn=dnorm(X*b);		
			z=X*b+(Yt-pn)%(1/dn);
			w=(dn%dn)/(pn%(1-pn));
			//cout << z << "m" << m;
			mat zst=rbind<mat>(z,m);
			//cout << "zst" << zst;
			mat sd=s.diag();
			mat wst=rbind<mat>(w,1/sd);
			mat W(n+p,n+p);
			W.fill(0);
			W.diag()=wst;	
			//cout << wst;
			//cout << W;
			mat r=inv(Xt.t()*W*Xt)*Xt.t()*W*zst;
			//cout << r;
			double temp=as_scalar(var(var(zst-Xt*r)));
			V=temp*inv(Xt.t()*W*Xt);
			//cout << "temp " << zst-Xt*r;
			//cout << V;
			mat foo=1/((r-m)-V.diag()-nu*s.diag());
			mat st=(1+nu)*(foo);
			//cout << st;
			//cout << "V" << V.diag();
			wst=rbind<mat>(w,1/st);
			W.diag()=wst;	
			r=inv(Xt.t()*W*Xt)*Xt.t()*W*zst;
			V=inv(Hessian(X,Y,r,s));
			e=abserr(r-b);
		//	cout << "\n"<< e << " ";
			b=r;
			i++;
			//cout << b;
			//cout << i;

		}	
		_m=b;
		_S=V;
	}
	inline mat pnorm(mat y)
	{
		int n=y.n_rows;
		mat x(n,1);
		for(int i=0;i<n;i++)
		{
			x(i,0)=cdf((*_d),y(i,0));
			if(x(i,0)==1)
			{
				x(i,0)=0.999999;
			}else if(x(i,0)==0){
				x(i,0)=0.000001;
			}
		}
		return x;
	}
	inline mat dnorm(mat y)
	{
		int n=y.n_rows;
		mat x(n,1);
		for(int i=0;i<n;i++)
		{
			x(i,0)=pdf((*_d),y(i,0));
		}
		return x;
	}


	inline double abserr(mat X)
	{
		double sum=0;
		int n=X.n_rows;
		for(int i=0;i<n;i++)
		{
			sum+=abs<double>(X(i,0));
		}
		return sum;
	}

	
	inline	mat Hessian(mat X,double *Y,mat theta,mat s)
	{
		int p=X.n_cols;
		mat sum(p,p);
		sum.fill(0);
		int n=X.n_rows;
		double xb=0;
		for(int i=0;i<n;i++)
		{
			
			xb=as_scalar(X(i,span::all)*theta);
			double pxb=pdf(*_d,xb);
			double gpxb=cdf(*_d,xb);
			mat xx=X(i,span::all).t()*X(i,span::all);	
			if(Y[i]==1)
			{	
				mat temp=(xb*pxb/gpxb+pow(pxb,2)/pow(gpxb,2))*xx;
				sum=sum+temp;				
			}else{
							
				mat temp=(-xb*pxb/(1-gpxb)+pow(pxb,2)/pow(1-gpxb,2))*xx;
				sum=sum+temp;				
			}
		}
		return sum+inv(s);

	}
	mat Get_Sig(void){return _S;}
	mat Get_Mu(void){return _m;}
	private:
	mat _m;
	mat _S;
	boost::math::normal_distribution<> *_d; 
};

